#include <CGalaxyComponent.h>




double CGalaxyComponent::getRadius(const CVector& rpos, CVector& P, double& dott)  {
  CVector p = rpos;
  dott = p.Dot(params->orientation);
  P = p - CVector(params->orientation*dott) ;
    if (abs(P.x*P.x)+ abs(P.y*P.y)+abs(P.z*P.z)<0.00001)
      return 0.00001;
  double rad = P.Length()/ params->axis.x;
  return rad;
}



double CGalaxyComponent::getHeightWeight(double height) {
  height  = 1.0-abs(height)*9.0;
  height = CMath::Minmax(height,0,1);
  height = pow(height, 2.0);

    return height;
}


double CGalaxyComponent::getPerlinCloudNoise(const CVector& p, double rad, double t, int N, double ks, double kadd) {
  CVector r = p;
  /*  if (rad<1.0)
      return 1.0;*/
  //  if (rad==1.0)
    {
      quaternion q;
      q = quaternion::from_axis_angle(params->orientation, abs(t)*CMath::sgn(params->orientation.y));
      RX.Identity();
      q.toMatrix(RX);
      r = RX.Mul(r);

     
  }
  double perlinnoise = 0;
  double max = 0;
  for (int ii=1;ii<N;ii++) {
    double k= ii*ks+kadd;
    double pp = abs(perlin->Get( (r.x + params->randShiftX*50.59)*k, (r.y+ 0.2 + params->randShiftY*50.12)*k, (r.z + params->randShiftY*49.87)*k))/k;
    //double pp = abs(perlin->Get( (r.x + 0.12 )*k, (r.y+ 0.2 )*k, (r.z + 0.74 )*k))/k;
    perlinnoise+=1.0/(pp+0.01);
    max += 6.0/ ((1.0/k) + 0.01);
    
  }
  average+=perlinnoise;
  count++;
  //cout << "avg:" <<average/count << endl;
  //cout << "max" << max << endl;
  return perlinnoise/max;//(average/count);
}


double CGalaxyComponent::findDifference(double t1, double t2) {
  double v1 = abs(t1-t2);
  double v2 = abs(t1-t2-6.28);
  double v3 = abs(t1-t2+6.28);
  double v4 = abs(t1-t2-6.28*2);
  double v5 = abs(t1-t2+6.28*2);

  double v = min(v1,v2);
  v = min(v,v3);
  v = min(v,v4);
  v = min(v,v5);

  return v;

}

double CGalaxyComponent::calculateArmValue(double rad, CVector& P, double armStrength) {


  /*    if (rad<0.001)
    return 1.0;
  */
  double v1 = getArm(rad, P, params->arm1);
  if (params->noArms==1)
    return v1;
  double v = max(v1, getArm(rad, P, params->arm2));
  if (params->noArms==2)
    return v;
  v = max(v, getArm(rad, P, params->arm3));
  if (params->noArms==3)
    return v;
  v = max(v, getArm(rad, P, params->arm4));

  return v;
}


double CGalaxyComponent::getArm(double rad, CVector p, double disp) {

  double t = getWinding(rad);
  double theta = -getTheta(p) + disp;
  double v = findDifference(t,theta);


  
  v = abs(v)/CMath::pi;

  //v*=pow((double)rad,0.1);  


  return pow(v,componentParams.arm);
}


double CGalaxyComponent::setupQuaternions() {
  // first, find the nearest plane
  CVector plane1(0,0,-1);
  CVector plane2(0,1,0);
  CVector plane3(1,0,0);
  double angle1 = params->orientation.Angle(plane1);
  double angle2 = params->orientation.Angle(plane2);
  double angle3 = params->orientation.Angle(plane3);
  // then choose the nearest one
  double angle = angle1;
  rA = &quatRot.x;
  rB = &quatRot.y;

  /*  cout <<"Orientation:" << params->orientation << endl;
  cout << "a1:" <<angle1 << endl;
  cout << "a2:" <<angle2 << endl;
  cout << "a3:" <<angle3 << endl;
  */
  CVector* plane = &plane1;
  if (angle2<angle) {
    plane = &plane2;
    angle = angle2;
    rA = &quatRot.x;
    rB = &quatRot.z;
  }
  if (angle3<angle) {
    plane = &plane3;
    angle = angle3;
    rA = &quatRot.y;
    rB = &quatRot.z;
  }
  //cout << "Plane chosen: " << *plane << endl;
  CVector dir = plane->Cross(params->orientation).Normalize();
  quat = quaternion::from_axis_angle(dir, angle);
  quat.toMatrix(rotmat);


  dir = plane1.Cross(params->orientation).Normalize();
  quat = quaternion::from_axis_angle(dir, angle1);
  quat.toMatrix(rotmatFull);

  // rotmatFull.Print();
  // After this, set up projection down to XZ-coordinates for twirl

}


double CGalaxyComponent::getTheta(const CVector& p) {
  //    cout << p << endl;

  quatRot = rotmat.Mul(p);
  return atan2(*rA, *rB);
}


double CGalaxyComponent::getWinding(const double& rad) {
  double r = rad;
  double t = atan((double)exp(-0.25/(0.5*r))/params->windingB)*2*params->windingN;
      
  double t2 = -1.0 * log(r*0.2);
	
  double scale = CMath::SmoothStep(0, 2*params->bulgeDust, r);

  scale = 0.2 + scale*0.8;

  if (params->innerTwirl == 0)
    scale = 1.0;

  //scale = 0.0;
  t = t*scale + t2*(1-scale);



  return t;

}



